function out = plotOutputFilter(outEstim)

setupFilter     = outEstim.setupFilter;
dateTimes       = outEstim.setupFilter.timeIndex;
timeIndexPlot   = round(year(dateTimes(1,:))/10)*10:10:round(year(dateTimes(end,:))/10)*10;

fontSizeValue   = 14-1;
figureSize      =  [0 0 1000+300 600+300]; % (0,0) size in x size in y

%termPrem        = outEstim.termPrem;
%% Model fit
figure('Name','Model fit: A prior','NumberTitle','off','Position',figureSize) 

dimy     = size(outEstim.outFilter.ybar,1);
dataUS   = setupFilter.data';
if dimy <= 9
    numCol = 3;
else
    numCol = 4;
end
numRow = ceil(dimy/numCol);
for i=1:dimy
    subplot(numRow,numCol,i)
    hold on
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),dataUS(i,setupFilter.numInitial+1:end),'-k')
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),outEstim.outFilter.ybar(i,setupFilter.numInitial+1:end),'-r')
    axis tight    
    datetick('x','yyyy','keeplimits','keepticks')   
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    title(setupFilter.dataLabels{i},'Interpreter','Latex','FontSize',fontSizeValue)
    if i == 1
        legend({'Data','Model'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'interpreter','Latex','FontSize',fontSizeValue);
    end
end

figure('Name','Model fit: Posterior','NumberTitle','off','Position',figureSize)
for i=1:dimy
    subplot(numRow,numCol,i)
    hold on
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
        dataUS(i,setupFilter.numInitial+1:end),'-k','LineWidth',1.0)
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
        outEstim.outFilter.yhat(i,setupFilter.numInitial+1:end)...
        -dataUS(i,setupFilter.numInitial+1:end)*0,'-r','LineWidth',1)
    plot(setupFilter.timeIndex(setupFilter.numInitial+1),-0.02)
    axis tight
    datetick('x','yyyy','keeplimits','keepticks')   
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    title(setupFilter.dataLabels{i},'Interpreter','Latex')
    title(setupFilter.dataLabels{i},'Interpreter','Latex','FontSize',fontSizeValue)
    if i == 1
        legend({'Data','Model fit'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'Position',[0.42732174840924 0.947971328952024 0.17616543107794 0.0413333309809365],...
               'interpreter','Latex','FontSize',fontSizeValue);
    end
end

figure('Name','Surveys: Posterior','NumberTitle','off','Position',figureSize)
dimy  = setupFilter.numMacro+length(setupFilter.maturities);
for i=1:length(outEstim.setupFilter.matConShortRate)
    subplot(2,2,i)
    hold on
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),dataUS(dimy+i,setupFilter.numInitial+1:end),'-k')
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),outEstim.outFilter.yhat(dimy+i,setupFilter.numInitial+1:end),'-r','MarkerSize',2)
    axis tight
    datetick('x','yyyy','keeplimits','keepticks')   
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    title(setupFilter.dataLabels{dimy+i},'Interpreter','Latex')
    title(setupFilter.dataLabels{dimy+i},'Interpreter','Latex','FontSize',fontSizeValue)
    if i == 1
        legend({'Data','Model'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'interpreter','Latex','FontSize',fontSizeValue);
    end
end

figure('Name','Surveys II: Posterior','NumberTitle','off','Position',figureSize)
indexSurvey = [1 4];
yMax        = 10;
yMin        = -0.1;
for i=1:length(indexSurvey)
    select = ~isnan(dataUS(dimy+indexSurvey(1,i),:));
    subplot(2,1,i)
    % NBER recession dates
    NBER        = outEstim.setupFilter.NBERrec(select);
    dates       = setupFilter.timeIndex(select);
    xaxis       = setupFilter.timeIndex(select);
    nn          = length(dates);
    idx         = nn:-1:1;
    reverseaxis = dates(idx);
    index       = 1:nn;
    reversIndex = nn:-1:1;
    grpyat = [xaxis       100*yMax*NBER(index,1);...
        reverseaxis NBER(reversIndex,1)*yMin*100];
    patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);    

    hold on
    plot(setupFilter.timeIndex(select),dataUS(dimy+index(1,i),select)*0,':k')
    tmp1 = plot(setupFilter.timeIndex(select),dataUS(dimy+index(1,i),select),'-k');
    tmp2 = plot(setupFilter.timeIndex(select),outEstim.outFilter.yhat(dimy+index(1,i),select),'-ok');
    axis([min(setupFilter.timeIndex(select)), max(setupFilter.timeIndex(select)),...
        -0.005 0.15])
    datetick('x','yyyy','keeplimits','keepticks')   
    set(gca,'FontSize',fontSizeValue+1)
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    title(setupFilter.dataLabels{dimy+indexSurvey(1,i)},'Interpreter','Latex','FontSize',fontSizeValue+3)
    if i == 1
        legend([tmp1(1), tmp2(1)],{'Data','Model'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'interpreter','Latex','FontSize',fontSizeValue+3);
    end
end


%% Analyzing the measurment errors
meanMerrors         = zeros(1,dimy);
stdMerrors          = zeros(1,dimy);
autoCorrMerrors     = zeros(10,dimy);
stdMerrorsNoOutlier = zeros(1,dimy);
numObsRemoveOutlier = nan(1,dimy);
for i=1:dimy
    meanMerrors(1,i)= nanmean(outEstim.outFilter.what(i,setupFilter.numInitial+1:end)); 
    stdMerrors(1,i) = nanstd(outEstim.outFilter.what(i,setupFilter.numInitial+1:end),0,2);
    tmp_i           = prctile(outEstim.outFilter.what(i,:),[2.5 97.5]);
    select          = (outEstim.outFilter.what(i,:) > tmp_i(1,1)).*(outEstim.outFilter.what(i,:) < tmp_i(1,2));
    select(1,1:setupFilter.numInitial) = 0;
    stdMerrorsNoOutlier(1,i) = nanstd(outEstim.outFilter.what(i,select==1));
    if i <= 9
        numObsRemoveOutlier(1,i) = length(select)-setupFilter.numInitial - sum(select);
    end
    select           = [zeros(1,setupFilter.numInitial)==1 ~isnan(outEstim.outFilter.what(i,setupFilter.numInitial+1:end))];
    if any(select > 0)
        tmp              = autocorr(outEstim.outFilter.what(i,select),10);
        autoCorrMerrors(:,i) = tmp(1,2:end)';
    end
end
nyMom      = outEstim.setupFilter.nyMom;
tmp        = outEstim.outFilter.what*outEstim.outFilter.what'/length(outEstim.outFilter.what);
corrMatrixMerrors = zeros(nyMom,nyMom);
for i=1:nyMom
    for j=1:nyMom
       corrMatrixMerrors(i,j) = tmp(i,j)/(sqrt(tmp(i,i))*sqrt(tmp(j,j)));
    end
end

figure('Name','Model fit: Measurement Errors','NumberTitle','off','Position',figureSize)
for i=1:dimy
    subplot(numRow,numCol,i)
    hold on
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
        outEstim.outFilter.what(i,setupFilter.numInitial+1:end)*10000,'-k')
    axis tight
    datetick('x','yyyy','keeplimits','keepticks')
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    %title([setupFilter.dataLabels{i},': std = ',num2str(round(stdMerrors(1,i)*10000)),'bp ',...
    %    'auto-corr(1) = ',num2str(round(autoCorrMerrors(1,i)*100)/100)],'Interpreter','Latex')
    %title([setupFilter.dataLabels{i},': std = ',num2str(round(stdMerrors(1,i)*10000)),' bps, $\rho_1 =$',' ',num2str(round(autoCorrMerrors(1,i),2))],'interpreter','Latex','FontSize',fontSizeValue)
    %title([setupFilter.dataLabels{i},': auto-corr(1) = ',num2str(round(autoCorrMerrors(1,i),2))],'Interpreter','Latex')
    title([setupFilter.dataLabels{i},': $\rho_1 =$',' ',num2str(round(autoCorrMerrors(1,i),2))],'interpreter','Latex','FontSize',fontSizeValue)
end

%% Analyzing the states and structural innovations
nx  = outEstim.model.nx;
nx1 = outEstim.model.nx1;
if outEstim.setupFilter.pruningOn == 1
    xhat          = outEstim.outFilter.xhat(1:nx,:);
    if setupFilter.appOrder == 1
        xhat(1:nx1,:) = xhat(1:nx1,:);
    elseif setupFilter.appOrder== 2
        xhat(1:nx1,:) = xhat(1:nx1,:) + outEstim.outFilter.xhat(nx+1:nx+nx1,:);
    elseif setupFilter.appOrder== 3
        xhat(1:nx1,:) = xhat(1:nx1,:) + outEstim.outFilter.xhat(nx+1:nx+nx1,:) + outEstim.outFilter.xhat(nx+nx1+1:nx+2*nx1,:);
    end
    xhatPruned    = outEstim.outFilter.xhat(:,:);
else
    xhat= outEstim.outFilter.xhat;
end
select = std(xhat')>1D-6;
[Axhat,~,Axhat_SE] = Var1(xhat(select,:)');
% Remove the intercept
Axhat = Axhat(1:end,1:end-1);
Axhat_SE = Axhat_SE(1:end,1:end-1);
Axhat_tstat = Axhat./Axhat_SE;

figure('Name','States: Posterior estimates','NumberTitle','off','Position',figureSize)
if nx <= 9
    numCol = 3;
else
    numCol = 4;
end
numRow = ceil(nx/numCol);
for i=1:nx
    subplot(numRow,numCol,i)
    hold on
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),xhat(i,setupFilter.numInitial+1:end),'-k')
    axis tight
    datetick('x','yyyy','keeplimits','keepticks')   
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    title(outEstim.model.labelx{i},'Interpreter','Latex','FontSize',fontSizeValue)
end

figure('Name','Structural Innovations: Posterior estimates','NumberTitle','off','Position',figureSize)
ne  = outEstim.model.ne;
nx  = outEstim.model.nx;
nx1 = outEstim.model.nx1;
T   = size(xhat,2);
inovShocks = zeros(nx,T);
for t=2:T
    if outEstim.model.pruningOn == 1
        xFitPruned_cu = stateTransition(xhatPruned(:,t-1),outEstim.model);
        xFit_cu = xFitPruned_cu(1:nx,1);
        if setupFilter.appOrder == 1
            xFit_cu(1:nx1,1) = xFit_cu(1:nx1,1);
        elseif setupFilter.appOrder == 2
            xFit_cu(1:nx1,1) = xFit_cu(1:nx1,1) + xFitPruned_cu(nx+1:nx+nx1,1);
        elseif setupFilter.appOrder == 3
            xFit_cu(1:nx1,1) = xFit_cu(1:nx1,1) + xFitPruned_cu(nx+1:nx+nx1,1) + xFitPruned_cu(nx+nx1+1:nx+2*nx1,1);
        end
    else
        xFit_cu = stateTransition(xhat(:,t-1),outEstim.model);
    end
    inovShocks(:,t) = xhat(1:nx,t)-xFit_cu(1:nx,1);
end
select_inov       = any(outEstim.model.eta>0,2);
inovShocksReduced = inovShocks(select_inov==1,:);
% Scaling inovShocks by their standard deviations
stdInov = outEstim.model.eta(outEstim.model.eta>0);
inovShocksReducedScaled = inovShocksReduced./repmat(stdInov,1,T);
meanInovShocksReducedScaled = nan(ne,1);
stdInovShocksReducedScaled  = nan(ne,1);
i = 0;
for j=1:ne
    if outEstim.model.eta(outEstim.model.nx1+j,j) > 0
        i = i + 1;
        tmp_i = prctile(inovShocksReduced(i,:),[2.5 97.5]);
        select = (inovShocksReduced(i,:) > tmp_i(1,1)).*(inovShocksReduced(i,:) < tmp_i(1,2));
        select(1,1:setupFilter.numInitial) = 0;
        meanInovShocksReducedScaled(j,1)  = nanmean(inovShocksReducedScaled(i,select==1),2);
        stdInovShocksReducedScaled(j,1)   = nanstd(inovShocksReducedScaled(i,select==1),0,2)...
            -0*mean(outEstim.outFilter.Smat(outEstim.model.nx1+j,outEstim.model.nx1+j,select==1));
    end
end
if ne <= 6
    numCol = 2;
else
    numCol = 3;
end
numRow = ceil(ne/numCol);
idx = 0;
for i=1:nx
    if select_inov(i,1) == 1
        idx = idx + 1;
        subplot(numRow,numCol,idx)
        hold on
        plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),inovShocks(i,setupFilter.numInitial+1:end),'-k')
        axis tight
        datetick('x','yyyy','keeplimits','keepticks')
        xticks(datenum(datetime(timeIndexPlot,1,1)))
        xticklabels(num2cell(timeIndexPlot))
        title(outEstim.model.labelx{i},'Interpreter','Latex','FontSize',fontSizeValue)
    end
end
select = std(inovShocks') > 1D-6;
[A,~,A_SE] = Var1(inovShocks(select,:)');
A_tstat = A./A_SE;
% Moments of inovShocks
meanInovShocks = mean(inovShocks,2);
stdInovShocks  = std(inovShocks,0,2);
corrInovShocks = corr(inovShocks');
nLagsInovShocks = 5;
autoCorrInovShocks = zeros(ne,nLagsInovShocks);
for i=1:ne
    tmp = autocorr(inovShocks(outEstim.model.nx1+i,:),nLagsInovShocks+1);
    autoCorrInovShocks(i,:) = tmp(1,2:nLagsInovShocks+1);
end



%% Plot of paistar
if 1 == 1
    figure('Name','Inflation shocks','NumberTitle','off','Position',figureSize)
    pospaistar      = find(strcmp(outEstim.model.labelx,'$\pi^*_{t}$'));
    paistar_cu      = outEstim.outFilter.xhat(pospaistar,:);
    pai             = setupFilter.data(:,4);
    paistar_cu      = mean(pai)+4*paistar_cu; %recentering around empirical mean
    corrM           = corr(pai,paistar_cu');
    %PHIpaiScaled_cu = (PHIpai_cu-repmat(nanmean(PHIpai_cu),length(PHIpai_cu),1))./repmat(nanstd(PHIpai_cu,0,1),length(PHIpai_cu),1);
    
    subplot(1,1,1)
    hold on
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),paistar_cu(setupFilter.numInitial+1:end),'-k')
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),pai(setupFilter.numInitial+1:end),'-r')
    hold off
    title(['$\pi^*_{t}$ and Inflation: corr = ', num2str(corrM)],'Interpreter','Latex')
    datetick('x','yyyy','keeplimits','keepticks')
    axis tight
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    set(gca,'FontSize',fontSizeValue)
    legend({'$\pi^*_{t}$ recentered','Inflation'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'interpreter','Latex','FontSize',fontSizeValue);    
   
end

%% Flot of real yields
% timeIndexPlot   = 1995:3:2020;  
% figure('Name','Model fit: Real Yields','NumberTitle','off','Position',figureSize) 
% for i=1:length(setupFilter.maturitiesReal)
%     subplot(2,2,i)
%     hold on
%     idx = setupFilter.maturitiesReal(1,i);
%     Tused = size(outEstim.outFilter.realYields(idx,1:end),2);
%     startIdx = max(find(isnan(setupFilter.realYields(1:end,i))))-4*3;
%     start2008 = min(find(year(setupFilter.timeIndex(1:Tused))>2007));
%     plot(setupFilter.timeIndex(startIdx:Tused),setupFilter.realYields(startIdx:Tused,i),'-k')
%     plot(setupFilter.timeIndex(startIdx:Tused),outEstim.outFilter.realYields(idx,startIdx:Tused),'-r')
%     axis tight
%     if ~isempty(start2008)
%         xline(setupFilter.timeIndex(start2008))
%     end
%     datetick('x','yyyy','keeplimits','keepticks')   
%     xticks(datenum(datetime(timeIndexPlot,1,1)))
%     xticklabels(num2cell(timeIndexPlot))
%     title(['$r^{(',num2str(idx),')}_{real,t}$'],'Interpreter','Latex','FontSize',fontSizeValue)
%     if i == 1
%         legend({'Data','Model'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
%                'interpreter','Latex','FontSize',fontSizeValue);
%     end
% end

%% Output
out.meanMerrors         = meanMerrors;
out.stdMerrors          = stdMerrors;
out.stdMerrorsNoOutlier = stdMerrorsNoOutlier;
out.numObsRemoveOutlier = numObsRemoveOutlier;
out.autoCorrMerrors     = autoCorrMerrors;
out.corrMatrixMerrors   = corrMatrixMerrors;
out.Axhat               = Axhat;
out.Axhat_tstat         = Axhat_tstat;
out.Ainov               = A(:,1:end);
out.Ainov_tstat         = A_tstat(:,1:end);
out.inovShocks          = inovShocks;
out.meanInovShocks      = meanInovShocks;
out.stdInovShocks       = stdInovShocks;
out.autoCorrInovShocks  = autoCorrInovShocks;
out.corrInovShocks      = corrInovShocks;
out.meanInovShocksScaled= meanInovShocksReducedScaled;
out.stdInovShocksScaled = stdInovShocksReducedScaled;
end

